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Abstract 


We developed a biochar model within the Agricultural Production Systems s[Mulator (APSIM) software that 
integrates biochar knowledge and enables simulation of biochar effects within cropping systems. The model has 
algorithms that mechanistically connect biochar to soil organic carbon (SOC), soil water, bulk density (BD), pH, 
cation exchange capacity, and organic and mineral nitrogen. Soil moisture (SW)-temperature—nitrogen limita- 
tions on the rate of biochar decomposition were included as well as biochar-induced priming effect on SOC min- 
eralization. The model has 10 parameters that capture the diversity of biochar types, 15 parameters that address 
biochar-soil interactions and 4 constants. The range of values and their sensitivity is reported. The biochar model 
was connected to APSIM’s maize and wheat crop models to investigate long-term (30 years) biochar effects on 
US maize and Australia wheat in various soils. Results from this sensitivity analysis showed that the effect of 
biochar was the largest in a sandy soil (Australian wheat) and the smallest in clay loam soil (US maize). On 
average across cropping systems and soils the order of sensitivity and the magnitude of the response of biochar 
to various soil-plant processes was (from high to low): SOC (11% to 86%) > N2O emissions (—10% to 43% 
43%) > plant available water content (0.6% to 12.9%) > BD (—6.5% to —1.7%) > pH (—0.8% to 6.3%) > net N min- 
eralization (—19% to 10%) > CO, emissions (—2.0% to 4.3%) > water filled pore space (—3.7% to 3.4%) > grain 
yield (—3.3% to 1.8%) > biomass (—1.6% to 1.4%). Our analysis showed that biochar has a larger impact on envi- 
ronmental outcomes rather than agricultural production. The mechanistic model has the potential to optimize 
biochar application strategies to enhance environmental and agronomic outcomes but more work is needed to 
fill knowledge gaps identified in this work. 
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Introduction 


Driven by the potential synergy of simultaneously 
sequestering carbon, enhancing soil quality and produc- 
ing renewable energy (Laird, 2008), research on the 
pyrolysis-biochar platform is growing rapidly and vari- 
ous initially small scale commercialization efforts are 
underway. Pyrolysis, the heating of biomass in the 
absence of oxygen yields bio-oil and syngas, which are 
bioenergy raw materials, and a solid biochar residue, 
which has potential use as a soil amendment (Laird 
et al., 2009). The physical and chemical properties of 
biochar, however, vary substantially depending on the 
nature of the biomass feedstock and the thermochemical 
processing conditions (Lehmann et al., 2011). Such vari- 
ation influences the stability of biochar carbon in soil 
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environments, hence carbon sequestration potential, as 
well as soil quality and agronomic responses to soil bio- 
char amendments. Furthermore, complex interactions 
between soil type, biochar type, climate, and manage- 
ment practices can substantially influence crop response 
to biochar applications. Thus a systems level under- 
standing of soil-crop-climate-biochar interactions is 
needed to predict a priori whether the application of a 
particular biochar on a particular soil will result in posi- 
tive, negative or neutral agronomic and environmental 
responses. 

To date, the majority of biochar studies have been car- 
ried out as short-term lab incubation or greenhouse pot 
experiments focusing on soil quality responses but 
neglecting crop and climate interactions. Results from 
these studies have demonstrated both desirable and 
nondesirable effects of various biochar types on soil 
properties and processes and plant growth. Biochar has 
been found to increase plant available water content 
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(PAWC), which is promising, given that soil water is 
often the main limiting factor of crop production (Basso 
et al., 2013). Most studies report a decrease in soil bulk 
density (BD) as a result of biochar application, which is 
another desirable characteristic as low BD enhances soil 
aeration and ultimately the ability of plant roots to 
explore and fully use the soil volume (Laird et al., 2010). 
The increase of soil organic carbon (SOC) as a result of 
biochar application is well documented in the literature 
(Brockhoff et al., 2010; Haefelea et al., 2011; Rogovska 
et al., 2014). However, biochar carbon is qualitatively 
different from biogenic SOC and thus is not expected to 
have the same influence on nutrient cycling and other 
soil processes as biogenic SOC. The effect of biochar 
additions on soil carbon and nutrient cycling is still 
under investigation and the experimental findings are 
inconclusive, as both increased and decreased green- 
house gas (GHG) emissions and N mineralization rates 
have been reported (Spokas & Reicosky, 2009; Clough 
et al., 2010; Bruun et al., 2011a; Strubel et al., 2011; 
Ameloot et al., 2013; Cayuela et al., 2014, 2015; Gai et al., 
2014; Case et al., 2015). Direct effects of biochar on soil 
cation exchange capacity (CEC), pH, PAWC, BD, and 
nutrient adsorption may have cascading influences on 
carbon and nutrient cycling. Furthermore, there is evi- 
dence that biochar amendments induce neutral, positive 
and negative priming effects on the turnover of SOC 
and the mineralization of crop residues (Cross & Sohi, 
2011; Zimmerman et al., 2011; Jones et al., 2012; Spokas 
et al., 2012; Wang et al., 2015). 

In contrast to incubation or pot studies, field trials 
examining the effect of biochar on intensively managed 
and fertile agricultural soils are currently rare and the 
effects on soil properties and processes are commonly 
less pronounced or not evident in these studies (Major 
et al., 2010; Haefelea et al., 2011; Guerena et al., 2013; 
Jones et al., 2012; Rogovska et al., 2014). Using a meta- 
analysis, Jeffery et al. (2011) found that biochar, on aver- 
age, increased crop yield by 15% in pot studies and by 
5% in field studies. The variation in yield response to 
biochar was high and included some negative 
responses, which reflects the complexity that exists in 
agroecosystems. 

The complexity in biochar research increases greatly 
when different types of biochar (feedstock source x pro- 
duction method x peak pyrolysis temperature) are con- 
sidered, which result in substantial variation in physical 
and chemical properties. Thus, determining where and 
what type of biochar to apply to optimize positive agro- 
nomic and environmental outcomes is very challenging. 
Yet this is essential for the strategic targeting of biochar 
applications and will add economic value to the pyroly- 
sis-biochar platform. To address this, biochar experi- 
mental findings must be integrated and re-analyzed 


under a common framework. Meta-analysis studies 
(Jeffery et al. 2011; Biederman & Harpole, 2013; Crane- 
Droesch et al., 2013; Cayuela et al., 2014; Liu et al., 2015; 
Wang et al., 2015) are the first steps in the integration of 
biochar information, but they lack the power of predic- 
tion that a mechanistic approach can offer. Our effort is 
to link the individual biochar effects in a way that incor- 
porates the complex interactions and feedback mecha- 
nisms. Jeffery et al. (2015) argued that biochar research 
needs a systems framework to allow effective compar- 
isons between experiments, across scales and locations, 
to guide policy and recommendations concerning bio- 
char application to soil. We hypothesize that agricul- 
tural systems models, which have a long history and 
have made significant contributions in agricultural 
research, can be used to quantify the complexities 
involved in the biochar-soil—plant—climate-management 
systems, and that doing so will help identify critical 
knowledge gaps, stimulate new strategic biochar 
research and accelerate commercialization. 

At present, modeling of biochar is still in its infancy, 
and published efforts have focused on a particular 
aspect of the system without addressing concurrent 
effects of biochar on soil water, nitrogen, carbon and 
their interactions (Woolf et al., 2010; Foereid et al., 2011; 
Woolf & Lehmann, 2012; Lychuk et al., 2015; Lim et al., 
2016). Our goal is to advance biochar modeling by 
accounting for concurrent biochar effects on soil carbon, 
water, nitrogen, and other soil properties and their 
dynamic interactions. Our specific objectives were: (i) 
develop a conceptual mechanistic biochar model by 
linking biochar effects to the soil-plant-atmosphere 
system within the Agricultural Production Systems 
sIMulator (APSIM) software platform (Holzworth et al., 
2014); (ii) provide a range of parameter values that 
address diversity among different biochar types and 
conduct a sensitivity analysis of those parameters; (iii) 
use the model to investigate long-term biochar effects on 
US maize and Australia wheat cropping systems; (iv) 
conduct a preliminary evaluation of the modeling 
approach by comparing predicted and measured soil 
and crop responses to biochar applications, and (v) iden- 
tify knowledge gaps and biochar research priorities. 


Materials and methods 


Biochar model development 


We synthesized biochar information and programmed 
the biochar model within APSIM (Holzworth et al., 
2014). APSIM is a worldwide leading agricultural sys- 
tems modeling platform, which is frequently updated, 
freely available, and contains more than 30 crop models 
and _ soil/environmental models along with flexible 
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specification of farm management. Hence it can assist 
investigation of biochar effects on different cropping 
systems growing in different environments (soil x cli- 
mate combinations). Because we developed and run the 
biochar model within APSIM we ensured that the bio- 
char model is compatible with the soil process models 
already developed within APSIM. 


Background on APSIM 


Agricultural Production Systems sIMulator is a field 
scale model that operates on a daily time step (www. 
apsim.info). Its modular design allows users to select a 
crop or a combination of crops in rotation, cultivars, 
soils, climates, and management practices. In turn the 
model provides outputs for numerous soil-plant vari- 
ables; including crop growth processes, soil water, soil 
temperature, nitrogen and carbon cycling, GHG emis- 
sions and residue dynamics (Probert et al., 1998; Thor- 
burn et al., 2001, 2010; Keating et al., 2003; Holzworth 
et al., 2014). The major processes simulated by the 
model and their feedback mechanisms are shown in 
Fig. 1. The model uses a multi-layer soil profile and soil 
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processes are simulated in each layer separately. SOC is 
divided into three pools, each with a different decompo- 
sition rate constant; including BIOM (microbial biomass 
pool with a fast decomposition rate), HUM (humic pool 
with an intermediate decomposition rate) and INERT 
(inert pool which does not decompose). Fresh organic 
matter (FOM) is a separate pool, which contains roots 
from the previous crop and crop residues if they have 
been incorporated into the soil by tillage. The FOM is 
also divided into three sub-pools reflecting fast 
(FPOOL1), slow (FPOOL2) and medium (FPOOL3) 
decomposition rates. Both FOM and BIOM are fast 
decomposing pools, but they have different carbon to 
nitrogen (C: N) ratios as they represent different 
organic materials. The C : N ratio of BIOM is fixed at 8, 
while the C : N ratio of FOM varies depending on the 
source of residue/roots. 


The biochar model 


Figure 1 provides an overview of the biochar model 
and how it integrates with nitrogen (N), and carbon (C) 
cycling in APSIM. The following sections describe the 
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Fig. 1 A simplified diagram illustrating APSIM’s soil carbon (top left panel), soil nitrogen (bottom left panel), soil water (top right 
panel) and crop (bottom right panel) modules. Boxes are state variable, solid arrows are rate variables and indicate material flow (e.g. 
carbon), broken arrows indicate information flow (e.g. priming), and circles are intermediate variables. Driving variables for the sys- 
tem (green circles at the top) include weather, soil, crop, and management. Shown in red are the state, rate, and intermediate variables 
that we hypothesize to be influenced by biochar amendments. RES, surface residue; FOM, fresh organic matter; BIOM, microbial pool; 
HUM, humic pool; INERT, inert pool; SAT, saturation point; DUL and LL, drained upper and lower limits; KL, parameter defining 
the ability of the roots to take up water; CEC, cation exchange capacity; and BD, soil bulk density. 
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underlined theory, the equations and the parameter 
values developed for the biochar model (Table 1). The 
development of each biochar parameter was accompa- 
nied by a sensitivity analysis. 


Biochar decomposition 


The starting point in the development of the biochar 
model was to define a new biochar C pool for APSIM 
that contains both labile and recalcitrant components as 
commonly described in the literature (Woolf et al., 
2010). The daily rate of biochar decomposition was 
modeled using a modified double exponential decay 
function (Archontoulis & Miguez, 2015) that accounts 
for water, temperature and N limitations: 


dltgc = BC, (1 = floss) 


fiabile ` Ino) 
365- MRT} ` 
- min(wr - tp- ne), 


(1 — fiabile) - Ino) 
365 - MRT2 


(1) 


` f carbon 


where, dltgc is the amount of biochar C that decom- 
poses every day (kg C ha’ day~'), BC, is the amount 
of biochar applied (kg dm ha), floss is the fraction of 
biochar that is lost during application (0 to 1), fearbon is 
the fraction of C in biochar (0 to 1), fiabile is the fraction 
of the labile biochar pool (0 to 1), MRT; and MRT; are 
the mean residence times for the labile and recalcitrant 
biochar pools, respectively, and wg ts and nę are SW, 
temperature and N modifiers (0 to 1) for the rate of bio- 
char decomposition. Given the lack of specific experi- 
mental data, we used modifiers similar to that used by 
APSIM for other organic materials including soil 
organic matter, residue and manure (Probert et al., 1998, 
2005; Thorburn et al., 2001, 2010; 2005). The equations 
and the parameters are provided below. 

Table 1 provides literature values for fcarbon, fiabile 
MRT; and MRT). Their values vary substantially among 
biochar types in relation to feedstock source and pro- 
duction method, and are therefore used in the model to 
characterize biochar diversity. These parameters are 
meaningful, measurable and user defined. The fearbon is 
determined by thermal combustion and may vary from 
0.1 to 0.9 (Spokas & Reicosky, 2009; Atkinson et al., 
2010; Brewer et al., 2011; Lehmann et al., 2011; Gai et al., 
2014). The fiabite is determined by proximate analysis 
{[(total C — fixed C)/total C] using ASTM standard 
method D 1762-84, 2007} and typically varies from 0.03 
to 0.30 (Woolf et al., 2010; Wang et al., 2015). In general, 
fiabile decreases with increasing pyrolysis temperature 
(Bruun et al., 2011b), but is also influenced by feedstock 
properties such as lignin content, and the thermochemi- 
cal process. The MRT; is based on incubation studies 
and varies from 0.3 to 25 years (Woolf et al., 2010; Wang 


et al., 2015). The MRT> can be estimated based on H : C 
elemental ratio determined by thermal combustion anal- 
ysis and varies from 50 to 50 000 years (Woolf et al., 
2010; Woolf & Lehmann, 2012). The parameter fioss is 
used to correct for losses that may occur during biochar 
application in the field. 

Equation (1) can take very different shapes depending 
on the parameter values used and the prevailing soil- 
weather conditions (Fig. 2). A sensitivity analysis of 
BC,, fearbony fiabile, MRT; and MRT, parameters showed 
that the MRT) has substantially less impact on the rate 
of biochar decomposition than the other parameters 
(Fig. 3a). The above parameters set the potential decom- 
position rate for biochar (Fig. 4). In practice, the prevail- 
ing soil-crop-weather conditions determine the actual 
decomposition rate, which is invariably lower. For 
example, Fig. 4 shows that the actual decomposition rate 
was 15% of the potential on an annual basis, with large 
day-to-day fluctuations. This is because soil water, tem- 
perature and N are rarely at optimum levels in the field; 
by contrast with laboratory incubations, which are com- 
monly run under optimum or near optimum conditions 
for biochar decomposition. As it is difficult to maintain 
optimum levels in the long term, this explains why the 
reported MRT; values in the literature declines as the 
experimentation duration increases (Wang et al., 2015). 

To account for limitations on the decomposition rate 
of biochar we used APSIM’s temperature and water 
modifiers on the rate of SOC decomposition (Probert 
et al., 1998). For the N modifier (n), we developed a 
new function similar to that currently used by APSIM 
for residue decomposition: 


—enrf_be - ((Soilgct /Navail) — Opt_bc) 
Opt_be , 


np = min (1 exp 
(2) 


where the cnrf_bc parameter determines the slope of 
the curve, Soilgcı is the labile biochar in the soil, Navai 
is the sum of inorganic N in the soil and N released 
during biochar decomposition (see below), and the 
opt_bc parameter determines the change from immobi- 
lization to mineralization. The default cnrf_bc and 
opt_bc values were assumed similar to APSIM residue 
(Table 1), but more research is needed to verify the val- 
ues of these constants. The sensitivity of the N modifier 
to these parameters is presented in Fig. 3c. The term 
Soilgct/Navai changes dynamically in relation to the 
amount of biochar in the soil and the soil mineral N 
available for biochar decomposition and in general this 
is highly variable (Fig. 4). The modifiers wz, tr and ny, 
though generic, are biochar specific because different 
biochar types affect the SW and N availability differ- 
ently. 


© 2015 The Authors. Global Change Biology Bioenergy Published by John Wiley & Sons Ltd., doi: 10.1111/gcbb.12314 


MODELING BIOCHAR EFFECTS ON SOILS AND CROPS 5 


Table 1 List of the parameters in the biochar model, with their definitions, units and the range of values found in the literature. 
BC1, BC2, BC3 are hypothetical biochar types used for sensitivity analysis and WC is a woodchip biochar used for model evaluation. 
Some important model variables are provided also 


Acronym Definition Unit Range BC1 BC2 BC3 wc 
Biochar type specific parameters and management 
BC, Biochar application rate Mg ha~! 0 to 150 20 20 20 0 to 98 
floss Biochar lost during application Otol Otol 0.02 0.02 0.02 0.02 
feitor Carbon fraction in biochar g g 0.17 to 0.88 0.7 0.6 0.8 0.78 
fiabile Fraction labile C in biochar Otol 0.03 to 0.3 0.2 0.05 0.3 0.13 
MRT, Labile pool mean residence time Years 0.1 to 30 3 10 1 1 
MRT, Recalcitrant pool mean residence time Years 50 to 50 000 200 150 2000 500 
BCgcEc Biochar effective cation exchange capacity cmolc kg~' 0 to 253 187 187 187 187 
BCLy Biochar liming value - 20 to 200 50 50 50 50 
CN_be Biochar C to N ratio = 7 to 769 100 50 300 132 
Depth Biochar incorporation depth mm 0 to 500 100 100 100 300 
Biochar and soil type specific parameters 
Prom Positive priming coefficient for FOM} m kg 'C Oto 0.138 0 0.1 0 0 
Pe Negative priming coefficient 1 for FOM m kg 'C  0to0.15 0 0 0.05 0 
Py Negative priming coefficient 2 for FOM m° kg“! C —0.30 to 0 0 0 —0.05 0 
Piom Positive priming coefficient for BIOM} m? kg! Ç 0 to 0.15 0 0.1 0 0 
Puum Positive priming coefficient for HUM§ m kg! C 0 to 0.15 0 0.1 0 0 
Pe Negative priming coefficient 1 for BIOM m° kg“! C —0.15 to 0 0 0 0.05 0 
Pe Negative priming coefficient 2 for BIOM m° kg C —0.15 to 0 0 0 —0.05 0 
Kads Langmuir adsorption coefficient mg ic 0 to 0.125 0.006 0.006 0.006 0.006 
Kaes Langmuir desorption coefficient mg L! 0 to 0.125 0.006 0.006 0.006 0.006 
Qu Quality modifier on LL Otol Otol 0 0 0 0 
Kput Slope parameter in DUL|| quality modifier Otol —1to0 —0.30 —0.15 —0.30 —0.15 
Kgp Slope parameter in BD** quality modifier sg —1to0 —0.15 —0.10 —0.15 —0.15 
BDa Fraction reduction in BD due to tillage Otol 0.6 to 1 0.8 0.8 0.8 0.8 
Upu Soil pH upper level = = 8.5 8.5 8.5 9.3 
Lou Soil pH upper level 7 = 3.5 3.5 3.5 3.5 
Biochar model constants 
ef bc Biochar carbon retention efficiency gg“ = 0.4 0.4 0.4 0.4 
fr_bcbiom Fraction of biochar to BIOM sg = 0.05 0.05 0.05 0.05 
cnrf bc Slope parameter in Eqn (2) s Otol 0.693 0.693 0.693 0.693 
opt_be CN ratio at or below biochar mineralizes N gg = 25 25 25 25 
Important model variables 
dltgc Daily rate of biochar C decomposition (kg C ha~! day~') 
dit_nbe Daily rate of N mineralized or immobilized during biochar decomposition (kg N ha7! day~') 
Soilgc Labile and recalcitrant biochar remaining in the soil (kg C ha~’) 
Soilgct Labile biochar remaining in the soil (kg C ha‘) 
We, te, Ne Water, temperature and nitrogen modifiers on the rate of biochar C decomposition (0 to 1) 


*Woodchip biochar used in Rogovska et al. (2014); see Fig 10. 
+FOM: fresh organic matter pool (mainly roots and residues). 


{BIOM: microbial soil organic matter pool (the labile portion of the soil organic matter). 


§HUM: humic soil organic matter pool (the bulk of the soil organic matter). 


ILL: soil water drained lower limit or permanent wilting point. 
||DUL: soil water drained upper limit or field capacity. 
*“BD: soil bulk density. 


Partitioning of the decomposed biochar carbon 


The gradual degradation of biochar (Eqn 1) emits some 
biochar C to the atmosphere, and transfers biochar C to 
the microbial biomass (BIOM) and humic (HUM) pools 
(Kramer et al., 2004; Brodowski et al., 2007; Kuzyakov 
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et al., 2009; Woolf & Lehmann, 2012). We modeled the 
partitioning of decomposed biochar C by defining an 


efficiency coefficient (ef_bc) that represents the fraction 
of decomposed biochar C that is retained in the soil and 
coefficients to distribute this retained C between the 
BIOM and HUM pools (Fig. 1): 
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Fig. 2 Simulated biochar decomposition over time as affected 
by the fraction of labile biochar (fiabile), and mean residence 
times for the labile and recalcitrant biochar pools (MRT; and 
MRT», respectively). Other biochar parameters are shown in 
the panel. The simulation refers to a maize crop growing in 
Ames, Iowa, under standard management practices for this 
region in a loamy soil (see Tables 2 and 3). 


dltgc = dlt_bc_CO; + dlt_bc_biom + dlt_bc_hum, (3) 


where: 


dlt_bc_CO; = dltgc - (1 — ef _bc), 
dlt_bc_biom = dltgc - ef_bc - fr_bcbiom, 
dlt_bc_hum = dltgc - ef_be.(1 — fr_bcbiom), 


where, dlt_bc_CO, is the daily loss of biochar C to the 
atmosphere (kg C ha”! day~'), dlt_bc_biom and dlt_be 
_hum are the daily flux of biochar C (kg C ha! day~') 
to BIOM and HUM, respectively, and fr_bcbiom is the 
fraction of retained biochar C that ends up in BIOM. 
The default value for ef_bc (Table 1) was assumed 
equal to the efficiency defined for humus decomposi- 
tion in APSIM (Probert etal., 1998) as has been 
assumed in other biochar models (Woolf & Lehmann, 
2012). The default value for fr_bcbiom was set so a 
small portion (0.05 x 0.4 = 2%; Table 1) of the decom- 
posed biochar C was partitioned to BIOM. This 
assumption is reasonable given experimental findings 
that the size of the microbial pool increases following 
biochar application (Clough et al., 2010; Bruun et al., 
2012; Zhang et al., 2014; Liu et al., 2015). The parameters 
described in this section have either not been experi- 
mentally measured or have been measured for only a 
limited range of biochar-soil-weather conditions (e.g. 
using isotopically labeled biochars, Zimmerman et al., 
2011). Figure 3b shows the sensitivity of the rate of CO) 
evolution to BC, (amount applied) and ef_bc (C reten- 
tion; Table 1). 


Biochar priming effects 


In addition to the changes in the size of the SOC pools 
(BIOM and HUM) imposed by biochar application and 
decomposition, some biochar types can alter the turn- 
over rates of the biogenic SOC pools (Fig. 1). This is 
known as the priming effect. Woolf & Lehmann (2012) 
investigated priming effects on the labile FOM, while 
Cross & Sohi (2011), Zimmerman et al. (2011) and Wang 
et al. (2015) studied priming effects on SOC. Using this 
information, we developed functions in the biochar 
model to address both positive (accelerates decomposi- 
tion rate) and negative (decelerates decomposition rates) 
priming effects. For the FOM pool, we defined a posi- 
tive priming coefficient (Prom) as a function of Soilge 
similar to Woolf & Lehmann (2012): 


Xpe =X. (1 + Prom - Soilgc /10 000), (4) 


where, Xp, is the new dynamic decomposition rate 
(day~') for each of the three APSIM FOM pools, X is 
APSIM’s default decomposition rate of each FOM pool 
(day—'), Prom is the positive priming coefficient 
(m? ke! C) that is assumed to be the same for all the 
FOM pools, Soilgc is the biochar remaining in the soil 
(kg C ha~’), and 10 000 is used for unit conversion. Bio- 
char negative priming effects on FOM were modeled by 
increasing the C transfer from FOM to the more stable 
HUM pool and simultaneously decreasing the fraction 
of FOM C that is respired; in other words the internal 
cycling of carbon was increased similarly to Woolf & 
Lehmann (2012): 


ef_fombc = ef_fom- (1+ P.-Soilgc/10000), (5) 


fr_fom_biombc = fr_fom_biom - (1 + P¢ - Soilgc/10 000), 
(6) 


where, ef_fombc is the modified ef_fom parameter that 
determines the amount of decomposed C retained in the 
system (0 to 1 value), fr_fom_biombc is the modified fr_- 
fom_biom parameter that determines the C flow from the 
FOM to BIOM pool (0 to 1) (note: the remaining C is allo- 
cated to HUM), and the parameters Pe and Py are the neg- 
ative priming coefficients (m? kg! C). Biochar induced 
positive and negative priming effects on SOC pools 
(BIOM and HUM) were computed in a similar way. 

All of the priming coefficients can be derived from 
lab incubation studies such as that of Kuzyakov et al. 
(2009) where isotopically labeled biochar was used to 
distinguish between CO, emission from biochar and 
SOC decomposition. Woolf & Lehmann (2012) reported 
that Prom (FOM positive priming; Table 1) varies from 
0% to 13.8% and the negative priming or stabilization 
coefficient varies from 0% to 3.6%. Estimates for Pgiom 
and Pyum (BIOM and HUM positive priming’s 
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Fig. 3 Sensitivity analysis of biochar model parameter. The output variables shown in each panel are: dltbc: daily rate of biochar C 
decomposition; dlt_bc_co2: daily rate of biochar C lost to the environment; dlt_nbe: daily rate of nitrogen mineralization /immobiliza- 
tion due to biochar decomposition; soil CEC BC: soil cation exchange capacity after biochar application; soil pH BC: soil pH after bio- 
char application; NH4 ads: soil NH4 adsorption due to biochar application; dlt_DUL: daily rate of change in drained upper limit; 
dlt_BD: daily rate of change in soil bulk density due to biochar application. The values in parentheses next to each parameter are the 
default values used in the sensitivity analysis. For symbols explanation and units see Table 1. 


coefficients) were difficult to extract from the literature alization (Cross & Sohi, 2011). Zimmerman ef al. (2011) 
because of the challenge in separating priming effects studied C mineralization rates for 20 biochar types on 
from the contribution of the labile biochar on C miner- five soils, and found that the magnitude of priming var- 
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Fig. 4 Simulated daily decomposition rate of biochar with their associated N immobilization (left panels) and their stress modifiers 
used to reduce the potential decomposition rate to actual (right panel; 1 means no stress; 0 means full stress). The biochar parameters 
values used are shown in Table 1 under bicohar type 1. The simulation refers to a maize crop growing in Ames, Iowa in a loamy soil 
(see Tables 2 and 3). The sharp change in the nitrogen stress factor on day of the year 125 is due to 180 kg N ha fertilizer application. 


ied from —52% to +89% at the end of the first year. 
Wang et al. (2015) in a review study reported lower bio- 
char priming values on SOC (from —8% to +20%) and 
that the direction and the magnitude of the priming 
depended on various factors including biochar charac- 
teristics, soil properties and experimental duration (bio- 
char age). The priming coefficients in our biochar 
model are linked to Soilgc (Eqns 4-6), which is 
determined by biochar characteristics (MRT, MRT, 
BC,, fearbory and fiabile) and soil properties (w,, tẹ and ng, 
see soil hydrology and BD later). In contrast, our model 
does not include yet an age factor because of the lack of 
specific data to support modeling. Nevertheless, poten- 
tial users can set different priming values in the model 
(Table 1) and test hypotheses related to the magnitude 
of the priming effect on soil processes and crop yields 
(see Results section). 


Biochar effects on soil N mineralization and 
immobilization 


During the decomposition of the biochar (Eqn 1), N is 
mineralized or immobilized depending primarily on the 
biochar C : N ratio, which is an input to the model 
(Table 1). Mineralization or immobilization of mineral 
N is determined as the balance between the release of N 
during biochar decomposition and the immobilization 
of mineral N during the formation of microbial biomass 
and new humic materials (Probert et al., 1998). The N 
needed for biochar decomposition is calculated as: 
dlt_bc_biom dlt_bc_hum 


dit_nbc_need = -N fiom? CN hum? O 


where, dlt_ nbc_need is the daily amount of N required 
for biochar decomposition (kg N ha! day™'), 


dlt_bc_biom and dlt_bc_hum is the daily biochar C that 
goes to BIOM and HUM pools, respectively (kg 
C ha! day), CN_biom is the C/N ratio of the BIOM 
pool (default = 8; apsim version 7.7), and CN_hum is the 
C/N ratio of the HUM pool (around 11-14; user 
defined). The N demand may be satisfied by the N that 
is released during the decomposition of the biochar, 
which is estimated as: 


dlt_nbc_released = dltgc/CN_be, (8) 


where, dlt_nbc_released is the mineralized N during 
biochar decomposition (kg N ha! day~') and CN_bc is 
the CN ratio of the biochar (Table 1; user defined). The 
C : N ratio of various biochar types varies from 7 to 769 
in the literature (Atkinson et al., 2010, Lehmann ef al., 
2011), and it is typically measured by thermal combus- 
tion. We believe that the C : N ratio of the labile fraction 
of biochar may be different from the C : N ratio of the 
recalcitrant fraction, but due to the lack of specific data 
and difficulties in measuring individual biochar pool C: 
N ratios, we currently use the C : N ratio of the whole 
biochar. This is another gap in biochar research. 

Finally, the net rate of N mineralization or immobiliza- 
tion during biochar decomposition was calculated as: 


dlt_nbc = dlt_nbc_release — dlt_bc_need, (9) 


where, dlt_nbc is the N mineralized or immobilized 
during biochar decomposition (kg N ha ! day‘). For 
example, in the case study shown in Fig. 4 approxi- 
mately 2 kg Nha! year! was immobilized during 
biochar decomposition. The mineralized or immobilized 
N is added or subtracted to/from the NH, pool on a 
daily basis, and this can greatly affect subsequent N 
cycling processes such as nitrification and denitrification 
as well as plant N uptake (Fig. 1). Our assumption here 
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is that a biochar with a CN ratio above 25 (parameter 
opt_bc; Table 1) will immobilize N, but the rate of N 
immobilization is generally very low because of the 
slow rate of biochar decomposition (Fig. 4). A sensitiv- 
ity analysis of the parameters affecting Eqn (9) showed 
that all are important (Fig. 3c) and can affect the process 
in very different directions. This varying sensitivity 
highlights the complexity that exists in the soil-biochar 
system. 


Biochar effects on soil CEC 


The biochar model has functions to estimate soil CEC 
with (Soilcgcpc) and without (Soilcgc) biochar present. 
When Soilcgc data are not available, the model esti- 
mates initial Soilcgc from SOC and noncarbonate clay% 
(naay) using a set of equations developed by Seybold 
et al. (2005) for different soil orders: 


Soilcec = f (naay, SOC, soil order), (10) 


The model then simulates the effect of biochar on soil 
CEC. Literature reports on the impact of biochar amend- 
ments on soil CEC are inconsistent with some reports 
showing an increase in CEC (Laird et al., 2010; Van Zwi- 
eten et al., 2010; Peng et al., 2011) and others showing no 
effect (Haefelea et al., 2011; Rogovska et al., 2014). Differ- 
ences are most likely related to the biochar type, the 
amount added, the native soil CEC, and biochar age. 
Over time, in soil environments surfaces of biochar parti- 
cles are oxidized forming various organic functional 
groups; dissolved organic compounds may be adsorbed 
on, or leached from, biochar surfaces; microorganisms 
may colonize biochar surfaces; carbonates and other ash 
components may be solubilized exposing new biochar 
surfaces; and the biochar equilibrates with the ambient 
soil pH. All of these aging processes influence the net 
impact of biochar amendments on soil CEC. We used a 
user defined effective biochar CEC (BCgcgc; Table 1), 
which is assessed as the net change in CEC per kg of soil 
approximately 1 year after a biochar addition relative to 
the mass of biochar in 1 kg of the soil. The BCgcgc 
approach integrates all of the above aging effects without 
identifying the specific mechanisms. The default value in 
the model for BCgcgc is 187 cmolc ke! (Laird et al., 
2010) while the range varies between 0.5 and 68 for fresh 
and 96 to 235 cmol kg”! for aged biochars, respectively 
(Cheng et al., 2008; Gai et al., 2014). The impact of a bio- 
char on soil CEC is calculated as: 


Soilcecge = Soilcec « (1 — Massg,) + BCecec - Masser, 
(11) 


where Soilcgc and Soilcgcgc are the CEC of the soil 
before and after a biochar addition, Mass,, is the mass 
fraction of biochar in the soil (g g estimated from 


layer depth, soil BD and biochar in the soil 
[Masse = Soilgc/(fearbon’ BDo ayer 10 000)I, and BCgEcEc 
is the effective CEC of the biochar (cmolc kg"). Sensi- 
tivity analysis showed that BCgcgc has a relatively 
small impact on Soilcrcge (Fig. 3d). 


Biochar effects on soil pH 


Most biochars are alkaline due to the presence of both 
carbonate rich inorganic ash, which is mixed with the 
biochar, and the presence of organic functional groups 
on biochar surfaces, both of which can react with pro- 
tons (Yuan et al., 2011; Fidel, 2012;). Hence the impact of 
a biochar addition on soil pH is similar to the addition 
of agricultural lime. We used a user defined lime value 
(BC, vy; Table 1) to represent the amount of acid that can 
be neutralized by the biochar. Typically, biochar lime 
values range from 20 to 200 cmol, ket. By contrast, 
pure CaCO; has a lime value of 2000 cmol, kg’. To 
estimate the impact of a biochar addition (or an agricul- 
tural lime addition) on soil pH we used a sigmoidal 
function (Nelson & Su, 2010) assuming that the buffer- 
ing capacity of soil is proportional to the soil CEC: 


; ; Massy, - BC 
SoilpHBc = Soilpy 4 Í as 


Soilcrcse | 
. (UpH = SoilpH) $ (SoilpH = Lpx) - Py (12) 
(UpH = Lpu) , 


where Soilpr and Soilpupc are the soil pH before and 
after biochar application, Mass¢, is the mass fraction of 
biochar in the soil (g py, BC,y is the biochar lime 
value (cmolc kg"), Upu and Lpy are the upper and 
lower limit of soil pH (here set at 8.3 and 3.5 assuming 
the agronomic pH range is limited by the precipitation 
of carbonates and solubility of aluminum hydroxide, 
respectively), and P; is a fitted constant that adjusts the 
buffering capacity of soils (here 10). Among these 
parameters, the Upu and the initial soil pH have the 
greater impact on Soilpypc (Fig. 3e). The soil pH 
changes rapidly after a biochar or lime addition or 
when the soil is amended with ammonium fertilizer 
otherwise the daily change in soil pH should be negligi- 
ble. Equation (12) can result in different shapes in 
response to input parameters (Table 1) and the obtained 
results are consistent with the magnitude of biochar 
effects on soil pH observed in the literature (Brockhoff 
et al., 2010; Clough et al., 2010; Major et al., 2010; Haefe- 
lea et al., 2011; Strubel et al., 2011; Yuan et al., 2011; 
Jones et al., 2012; Jien & Wange, 2013). 


Biochar effects on NH; adsorption and desorption 


The adsorption of the NH, due to biochar addition is 
considered to be mainly due to the CEC. We quantified 
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the effect of biochar on NH, adsorption using the 
Langmuir model, which is frequently used for describ- 
ing sorption isotherms (Gai et al., 2014): 


Kaas « (Soilcecgc /Soilcsc) 


NHas = NH3- - - 
sads 4 1+ Kaas - (Soilcecse / Soilcec) 


> (13) 


where NH4zaas is the adsorbed NH; from the soil solu- 
tion (mg L~+), NH; is the concentration of the ammo- 
nium in the soil solution (mg L), Kaas refers to 
adsorption capacity (L mg ') and is defined by the user 
in the biochar model (Table 1; Fig. 5), Soilcrcpc and 
Soilcgc is the soil CEC before and after biochar applica- 
tion. Literature values for Kaqs vary from 0.0026 to 
0.125 L mg | (Table 1; Ding et al., 2010; Gai et al., 2014). 
Equation (13) is very sensitive to the Kaqs parameter 
(Fig. 3f). Equation (13) shows that the ability of the soil 
to adsorb NH, due to the biochar addition is propor- 
tional to CEC, which is a function of the added biochar 
in the model (see Eqn 11). 

A CEC of 1 cmol ke! is equivalent to a maximum 
sorption of 0.14 mg NH4 g! biochar, which means that 
biochar’s maximum adsorption capacity is from 0.07 to 
9.5 mg NH; gp" biochar assuming fresh biochar CEC 
values from 0.5 to 68 (Gai et al., 2014; Jassal et al., 2015). 
However, recent literature findings have shown even 
higher maximum adsorption capacities, i.e. 1.1-28.4 (Jas- 
sal et al., 2015) and 2.6-15.4 (Gai et al., 2014), which led 
to a hypothesis that most of the sorption must be due to 
physical entrapment of NH; in (fresh) biochar pores 
(Spokas et al., 2011; Jassal et al., 2015). The above values 
clearly show that the ability of the biochar to adsorb 
NH; is not universal, but depends on several factors 
including biochar type, age, and soil type (Ding et al., 
2010; Nelissen et al., 2012; Yao et al., 2012; Gai et al., 
2014; Sarkhot et al., 2014; Jassal et al., 2015). Equa- 
tion (13) can be extended further to account for more 
mechanisms or variables such as O%, C%, O/C ratio, 
and H/C ratio, which have been found to correlate well 
with the rate of NH, adsorption (Gai et al., 2014). 

In contrast to NH, adsorption, much less information 
is available about the release of the adsorbed NH; to the 
soil solution. If the soil-biochar system adsorbs NH,, 
then theoretically, at some point in time it should release 
it to the soil solution until equilibrium is achieved. Based 
on this assumption, we modeled NH, desorption simi- 
larly to adsorption (see Eqn 13). The rate of desorption is 
driven by a user defined Kaes parameter (Table 1; Fig. 5). 
Figure 5 illustrates Kaas and Kaes effects on simulated 
NH; over time. The value of Kaes is largely unknown 
and we used in the biochar model similar values as the 
Kaas, despite the fact that some short duration experi- 
ments have shown lower values. Clearly more research 
is needed in this area given the importance of NH; in N 
cycling (Sarkhot et al., 2014; Jassal et al., 2015). 
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Fig. 5 Simulated soil NH,-N (mg L“; 0-10 cm depth) 
response to adsorption (Kaas) and desorption (Kaes) coefficients. 
Biochar type 1 parameter values were used in this simulation 
(see Table 1) expect for the Kaas and Kaes parameters, which val- 
ues are shown in the panel. The simulation refers to a maize 
crop growing in Ames, Iowa in a loamy soil (see Tables 2 and 
3). The sharp increase in NH4-N on day of the year 123 is due to 
application of 100 kg N ha7! urea-N fertilizer in this example. 


Biochar effects on soil water and BD 


The effect of biochar on APSIM soil water retention 
parameters [lower limit (LL, mm mm _ '), drained upper 
limit (DUL, mm mm’), and saturation point (SAT, 
mm mm_})J, as well as soil BD (in g cm™°) was modelled 
by incorporating relationships between SOC and these 
parameters from the Saxton & Rawls (2006) pedotransfer 
functions into the biochar model (Fig. 1). Inputs for the 
model are clay and sand content of the soil and the rate of 
change in SOC due to biochar addition. The rate of change 
in SOC can be positive as a result of biochar addition or 
negative when SOC declines over time. Clay and sand 
contents mediate the rate of change. Outputs from the 
model are the rate of change in LL, DUL and BD on a 
daily basis. The rate of change in BD is used to compute 
SAT from the change in total porosity which is a function 
of BD [SAT = (1 — BD/2.65) x 0.95]. The value of 2.65 is 
the particle density of the mineral soils (g cm” °) and 0.95 
reflects the air filled porosity (Dalgliesh & Foale, 2005). 
The Saxton & Rawls (2006) equations are driven by 
the amount of SOC and have not been tested for 


© 2015 The Authors. Global Change Biology Bioenergy Published by John Wiley & Sons Ltd., doi: 10.1111/gcbb.12314 


MODELING BIOCHAR EFFECTS ON SOILS AND CROPS 11 


biochar. We hypothesized that the effect of biochar on 
soil hydrology will be somewhat different from that of 
SOC as biochar C is qualitatively different from SOC. 
Literature supports this hypothesis as it shows that dif- 
ferent types of biochars have very different surface pore 
spaces (0.2 to 8.2 nm; Gai et al., 2014; Lim et al., 2016) 
and variable effects on DUL, LL, and BD, which cannot 
be explained by the amount of biochar applied alone 
(Brockhoff et al., 2010; Laird et al., 2010; Liu et al., 2012; 
Basso et al., 2013; Cornelissen et al., 2013; Herath et al., 
2013). These studies showed that the LL is rarely chan- 
ged due to biochar addition, DUL and SAT can be chan- 
ged but the rate of change depends on the initial SOC 
before application and on biochar type. To account for 
the literature findings and to add further flexibility in 
the model, we developed and incorporated into the bio- 
char model 0-1 quality modifiers on the rate of change 
in LL, DUL and BD (Qır, Qput, and Qep, respectively). 
The Qı; was set to zero (Table 1) because biochar does 
not change LL, in contrast to Saxton & Rawls (2006) 
equations that show a change in LL. The Qpur was set 
as a function of SOC [Qpur = exp (—KpuL:SOC)]. The 
Kpu parameter is used to correct for the negligible 
effect of biochar on DUL in soils rich with organic C. 
The Qgp was set as function of SOC similar to Qpur but 
with a different input parameter (see Kgp in Table 1) 
because the effect of biochar on SAT is usually larger 
than that on DUL (e.g. Liu et al., 2012). The parameters 
Qır, Kpur and Kep are user defined (Table 1). By setting 
Qır to 1 and Kpur and Kgp to zero, the quality modi- 
fiers are cancelled, and the equations take the original 
Saxton & Rawls (2006) form. A sensitivity analysis of 
the biochar hydrological parameters is presented in 
Fig. 3g. 

The biochar model is capable of representing the 
effect of soil type (via texture), the physical aspects of 
biochar (amount added) as well as quality aspects of 
biochar (pore space) on soil water retention and BD. As 
SOC increases due to biochar addition, DUL and conse- 
quently PAWC increase and the magnitude is larger in 
coarse textured than in fine textured soils. Figure 6 
shows simulated PAWC response to biochar additions 
in soils with different levels of organic C. The magni- 
tude of the simulated changes in PAWC is consistent 
with literature findings (Brockhoff et al., 2010; Laird 
et al., 2010; Liu et al., 2012; Basso et al., 2013; Cornelissen 
et al., 2013; Herath et al., 2013). 

For the biochar effect on BD we considered a tillage 
factor in addition to the effect of C. This is because most 
laboratory incubation studies have shown a decreased 
BD in biochar amended soils relative to controls (Laird 
et al., 2010; Lim et al., 2016), but results from field stud- 
ies are more equivocal with most studies either showing 
no significant effect or a decrease in BD after biochar 


application (Haefelea et al., 2011; Rogovska et al., 2014). 
Under field conditions biochar is commonly incorpo- 
rated by tillage and it can be difficult to distinguish the 
biochar effect from the tillage effect on soil BD. To 
address this challenge we used the law of maximum 
absolute change between the two processes: biochar 
effect on BD (through the Saxton and Rawls pedotrans- 
fer functions) and tillage effect on BD (through Eqn 14). 
To simulate tillage effects on BD we used the following 
equation (Andales et al., 2000): 


BD = BD, + (BD, — BD,) - ASO”), (14) 


where BD is the calculated BD (g cm~°) that changes 
over time. BD, is the user input initial BD (g cm~’). 
BD, is the BD just after tillage operation (this value is 
defined by the user and varies with tillage equipment; 
Table 1). The AS parameter reflects soil aggregate sta- 
bility and equals 5-(1-0.205 - SOC), and Crain is the 
cumulative amount of precipitation since the latest til- 
lage event (see Andales et al., 2000). Figure 3h illus- 
trates the sensitivity of BD related parameters on 
Eqn (14). In general, the biochar model decreases BD 
after tillage to BD, (Fig. 7; Table 1). Then, with precipi- 
tation as the main driver, BD returns to the initial con- 
dition or to the new condition set by the added biochar 
C (Fig. 7). Literature findings have shown an average 
reduction in soil BD of 11.1 + 8.6% due to biochar 
addition (Lim et al., 2016). The reduction in BD due to 
biochar in this example was 5.8% (Fig. 7; loamy soil). 
The final BD estimate from the biochar model is used 
to compute SAT on a daily basis. 

It should be mentioned that the dynamic changes in 
LL, DUL, and SAT have a substantial impact on the sim- 
ulation process because they are linked to water modi- 
fiers that can alter C and N cycling and plant growth 


10 T-~o-SOC (0.5) 
-CESOC (1.0) 
T-soc (2.0) 
6 | ~-SOC (3.5) 


Simulated change in PAWC (%) 


0 20 40 60 80 
Biochar application rate (Mg ha’) 


Fig. 6 Simulated changes in plant available water content 
(PAWC = DUL - LL) by the biochar model as affected by the 
application rate and soil organic carbon. A biochar type 1 
parameter values were used in this simulation (see values in 
Table 1; expect the BCo which varied) and a loamy soil (see 
Table 2, expect the initial soil organic carbon which varied). 
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(Figs 1 and 4). As the biochar model affects both DUL 
and SAT, the SW will change as well. This in turn affects 
the water filled pore space (WFPS = SW/SAT), which is 
one of the main drivers for N20 emissions in the model. 
This shows the feedback mechanisms that exist in the 
model (Fig. 1). It would be possible to incorporate bio- 
char effects on saturated hydraulic conductivity (Ksat) in 
the model through the Saxton and Rawls equations, but 
it was not finally included in this version for two 
reasons: (i) literature findings are still contradictory (see 
Table 3 in Lim et al., 2016) and (ii) a preliminary sensitiv- 
ity analysis of Ksat (i.e. 100 mm day ' + 50% change) 
had a very minor impact on the simulation output (data 
not shown). 


Summary of biochar parameters 


We classified the biochar model parameters in three 
major categories (Table 1): (i) those that are highly 
dependent on biochar type (production 
method x feedstock source x production temperature) 
and management (application amount, incorporation 
depth); (ii) those that are determined by soil and bio- 
char interactions (effects on soil hydrology, priming); 
and (iii) finally those that are assumed to be constants 
in the model. All of the above parameters are impor- 
tant and affect various soil processes and properties 
differently (Figs 3 and 8). Parameters for the first cate- 
gory can be extracted from various literature 
resources or measured as outlined above or can be 
taken from the Biochar Engineering web tool (http:// 
spark.rstudio.com/veromora/BiocharEng/; Morales 
et al, 2015) or other biochar databases (http:// 
biochar.ucdavis.edu/). For the second and third cate- 
gory, literature information is incomplete and specific 
experiments are needed to fill modeling knowledge 
gaps (see Discussion later). 


Preliminary model testing 


We used experimental data from Rogovska et al. (2014) 
to test various aspects of the biochar model. Briefly, in 
this experiment a woodchip biochar (see parameter 
values in Table 1) was applied at six rates (range from 0 
to 98 Mg ha‘) to a maize field located in central Iowa. 
For more details we refer to Rogovska et al. (2014). Our 
simulation was for a continuous maize crop (Table 2) 
growing in a Clarion soil (Table 3) in Ames, IA. 


Results 


Long term sensitivity analysis on systems performance 


A sensitivity analysis of the individual biochar model 
parameters on key soil processes was presented above 
(Fig. 3). However, this type of sensitivity analysis (Fig. 3, 
change one parameter at a time while keeping all others 
constants) has some important limitations when the 
results are to be scaled up to the final product (i.e. grain 
yield; Fig. 1). This is because the results of the sensitivity 
analysis are dependent on the soil type used, cropping 
system, climate, management, and biochar parameters. 
The fact that there is not yet a biochar type that can serve 
as a reference in the sensitivity analysis, makes the initial 
choice of the biochar parameters very critical. 

Given all these issues we performed a second 
sensitivity analysis in which we explored the effect of dif- 
ferent biochars (named BC1, BC2, BC3; Table 1) on 
diverse cropping systems (crop, management, climate 
and soils; Tables 2 and 3). In this analysis we output sev- 
eral production and environmental variables in the long- 
term (30 years) to get a complete picture of systems per- 
formance. The results are presented in Fig. 8. On aver- 
age, across five cropping systems and three biochar 
types, the order of sensitivity of various soil-plant 
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Fig. 7 Simulated effects of biochar and tillage on soil bulk density (left panel) and cumulative rainfall as tillage (right panel). The 


biochar applied on day of the year 96 and incorporated by tillage on the same date. The biochar parameter values are shown in 
Table 1 (biochat type 1). The simulations refer to a maize crop growing in Ames, Iowa in a loamy soil (see Tables 2 and 3). 
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Fig. 8 Simulated long term biochar effects on soil-crop-atmospheric variables of two cropping systems (US maize and AU wheat) 
growing in different soils and climatic conditions (see Tables 2 and 3) as affected by three different biochars (see Table 1 for the 
parameter values). The duration of each simulation was 30 years and the results averaged across this period (expect the value of 
SOC) and compared to the control simulation (without biochar) for each cropping system. A negative value indicates that the biochar 
reduced the output variable (yield, CO, emissions, etc.) and vice versa. The values in each panel and cropping system indicate the 
long term value of the control simulation. 


Table 2 Characterization of the soils used in the simulation experiment. Nicollet, Clarion and Zenor soils are found in Midwest US, 
while Red clay and yellow sand are typical in Western Australia 


Soil profile Top soil layer where biochar has been incorporated 

Soil Total DUL LL 

depth PAWC OC (mm (mm Clay Sand BD pH 
Soil name Soil type (m) (mm) (%) mm!) mm!) (%) (%) (g cm) (-) 
Clarion Loam 1.82 248 2.03 0.292 0.149 21.0 41.6 1.37 6.5 
Nicollet Clay loam 1.54 218 3.19 0.305 0.174 27.5 38.1 1.20 5.5 
Zenor Sandy loam 1.44 118 1.63 0.182 0.096 12.5 67.9 1.52 6.4 
Red clay Clay 1.80 140 1.11 0.370 0.180 40.0 20.0 1.30 8.8 
Yellow sand Loam sandy 2.50 111 0.48 0.105 0.045 10.0 90.0 1.62 5.7 


PAWC, plant available water content; DUL, drained upper limit or field capacity; LL, drained lower limit or permanent witling point; 
BD, bulk density. 


processes and properties to biochar application was ization > CO, emissions > WFPS > grain yield > crop 
(from high to low): SOC > N20 emis- biomass. The direction and the magnitude of the change 
sions > PAWC > soil BD > soil pH > soil net N mineral- (compared to the control) was highly variable among 
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Table 3 Characterization of climate and management practices used in the simulation experiment 


Region and cropping system 


Midwestern US maize 


Western Australia wheat 


Site name Ames, Iowa Wubin 

Coordinates 42.0°N, —93.8°W —30.1°N, 116.6°W 

Frost free days (T min > 0 °C) 223 365 

Mean annual temperature (°C) 9.3 19.1 

Mean annual rainfall (mm) 900 315 

Crop Maize Wheat 

Cultivar name B_110* Mace 

Planting date May 5th every year Variable between Apr 25th to Jul 15th 
Plant density (pl m~’) 9 100 

Row spacing (cm) 70 30 

Sowing depth (cm) 5 3 

Tillage Every year Once during biochar application 
Tillage type Conventional Conventional 

N fertilization date At planting At planting and 40 days after 

N rate (kg N ha7') 180 20 + 40 = 60 

N type Urea_N Urea_N 


*It refers to a 110 day hybrid (see specifics in Archontoulis et al., 2014). 


cropping systems and biochar types: SOC varied from 
11% to 86%, N20 emissions from —10% to 43%, PAWC 
from 0.6% to 12.9%, soil BD from —6.5% to —1.7%, soil 
pH from —0.8% to 6.3%, soil net N mineralization from 
—19% to 10%, CO» emissions from —2.0% to 4.3%, WFPS 
from —3.7% to 3.4%, grain yield from —3.3% to 1.8%, and 
crop biomass from —1.6% to 1.4%. The effect of biochar 
on various variables was largest in the loam sandy soil 
(Australian wheat) and smallest in clay loam soil (US 
maize) (Table 2, Fig. 8). 

Among the tested variables, SOC, PAWC and BD 
showed a monotonic increase or decrease while all the 
other variables showed both increased and decreased 
trends over time (Figs 8 and 9). Our results indicated 
that the effect of biochar on crop production variables 
was much smaller compared to the effect on environ- 
mental variables such NO and CO, emissions. The 
reason might be the feedback mechanisms that exist in 
biochar-soil-crop-atmosphere systems and are inte- 
grated into our model. These mechanisms may sub- 
stantially counterbalance biochar effects on individual 
processes or the biochar effects are not large enough to 
affect crop yields. For example, an increase in SW that 
favors crop growth might trigger NO3-N loss via deni- 
trification or leaching, which can limit crop growth. 
Given that crop yield is the final product of several 
interrelated soil-plant processes (Fig. 1), our model 
analysis suggests that the effect of biochar diminishes 
all the way to the final product, grain yield. This result 
also shows that it is possible for biochar to decrease 
GHG emissions without negatively impacting crop 
yields which is important (see the magnitude of the 
biochar response: —10% to 43% for NO vs. —3.3% to 


1.8% for grain yield). This can be achieved by identify- 
ing appropriate biochar types within a cropping system 
and environment. The biochar model coupled with 
appropriate optimization functions can assist in this 
task. 

The magnitude of the obtained results in Fig. 8 is con- 
sistent with literature findings on biochar and reveal for 
the first time the trade-offs that exist in the biochar— 
soil-crop—atmosphere system. For example, the applica- 
tion of BC2 resulted in higher grain yield and biomass 
production and at the same time in higher GHG emis- 
sions. That was caused by the higher PAWC and soil N 
availability (see PAWC and N mineralization in Fig. 8). 
The decrease in WFPS was not enough to reduce N?O 
emissions. In contrast, the BC3 produced the opposite 
responses compared to BC2. It decreased both grain 
yield and GHG emissions variables. Literature studies 
have shown decreased and increased NO and CO, 
emissions due to biochar application (Clough et al., 
2010; Bruun et al., 2011a; Case et al., 2012, 2014, 2015; 
Ameloot et al., 2013; Cayuela et al., 2014, 2015; Liu et al., 
2015). Our study adds significant value to the biochar 
research because it reveals for the first time the reasons 
responsible for increased or decreased GHG emissions 
and grain yields (see parameters in Table 1, mecha- 
nisms in Fig. 1 and the magnitude of the response in 
Fig. 8). 

Figure 8 presents the average effect of the biochars 
across 30 years of cropping. It should be mentioned that 
the response of biochar on production and environmen- 
tal variables was not consistent across years, but vari- 
able with positive, neutral and negative responses over 
time (Fig. 9). This is consistent with observations and 


© 2015 The Authors. Global Change Biology Bioenergy Published by John Wiley & Sons Ltd., doi: 10.1111/gcbb.12314 


MODELING BIOCHAR EFFECTS ON SOILS AND CROPS 15 


107 _o Grain yield 


L == CO2 emissions 


- N20 emissions 


% change from control (no biochar) 
e 


6 8 10 12 14 16 18 
Years after biochar application 


Fig. 9 Simulated temporal variability of biochar response 
(BC1; Table 1) to grain yield, annual cumulative CO, emission 
and annual cumulative N20 emissions as compared to the con- 
trol (no biochar). The simulation refers to US maize growing 
on a loamy soil (see Tables 2 and 3). A negative value indicates 
that the biochar reduced the output variable (yield, CO, emis- 
sions, etc.) and vice versa. 


highlights the importance of climate and the carry-over 
effects on the biochar response. This also verifies the 
need for long term field assessments of biochar before it 
can be recommended for large-scale applications (Jones 
et al., 2012; Nelissen et al., 2012). Modeling can assist in 
this task by extending knowledge, but long-term field 
trials are ultimately needed to calibrate and validate 
modeling. 


Preliminary biochar model evaluation 


Figure 10 shows the goodness of fit of the biochar 
model to several soil-plant variables obtained under 
field conditions (Rogovska et al., 2014). The calibrated 
biochar parameter values used in this simulation are 
shown in Table 1. Overall, the model matched experi- 
mental observations with a mean relative absolute error 
from —0.4% to 13.1%, which is very good. The predic- 
tion of soil BD and pH was the most accurate, while the 
ability of the model to predict more complex variables 
such as SOC was lower (Fig. 10). In this particular set of 
data (year: 2012, location: Ames, IA), the model pre- 
dicted that corn grain yields were limited by SW, which 
agrees with the experimental observations (Rogovska 
et al., 2014). The model predicted that biochar applica- 
tion positively influenced SW (Fig. 10f), but this 
increase was not enough to compensate for the extre- 
mely dry conditions that occurred in 2012. 


Discussion 


The potential of the biochar model 


In a recent state-of-the-art biochar review study, Jeffery 
et al. (2015) argued for the need of a mechanistic under- 
standing of the effects of biochar application in order to 


quantify and predict benefits and trade-offs. In this 
paper we presented the first complete biochar simula- 
tion model that integrates and synthesizes broad bio- 
char knowledge and enables simulation of biochar 
effects on diverse cropping systems growing in different 
environments. We also demonstrated the potential of 
the model to assist system level understanding and 
assessments of biochar and their trade-offs (Figs 8 and 
10). The biochar model can assist biochar research in 
various ways. It can be applied to analyze data from the 
vast lab incubation studies or field experiments and to 
extrapolate the results beyond the study period in time 
and space in order to evaluate long-term consequences, 
which are definitely needed. Once the model fits the 
experimental observations (calibration task), it can pro- 
vide information on soil-plant processes that have not 
been measured experimentally and therefore allow for a 
system analysis of the biochar effects. The model can be 
used to test or generate hypotheses to inform experi- 
mental research or to perform experiments that are 
expensive, time consuming and difficult to implement 
in real-life situations. For example, what is the optimum 
biochar rate and type to maximize yields and minimize 
GHG emissions from a specific cropping system in a 
particular location? The model can provide an answer 
to this question relatively quickly and inexpensively. By 
running a factorial simulation with multiple biochar 
rates, multiple biochar types within a rate, different 
management practices (N application rate, cultivars, 
planting dates, tillage type and time) for each biochar 
type, and considering 30 years of weather variability, 
the system can be optimized and the win-win combina- 
tion can be found theoretically. Then the best scenario 
can be tested experimentally to verify model predic- 
tions. The biochar model can be coupled with other 
advanced simulation platforms such as pSIMS (Elliot 
et al., 2014) to assist regional- to global-scale assess- 
ments of biochar. 


Current status and future steps 


As this was the first attempt to develop a generic tool 
for system assessments, we initially focused on develop- 
ing algorithms and predicting the direction and the 
magnitude of the biochar responses (Figs 4-8). Our 
work on model conceptualization and development is 
driven by diverse literature information and not from a 
site-crop-biochar specific experimental dataset. Usually 
modelers build their models using a particular experi- 
mental dataset and use another data set for validation 
(Lychuk et al., 2015). In contrast, we developed a com- 
prehensive theoretical framework by integrating a wide 
range of biochar research findings and synthesized it 
from a systems perspective and performed a 
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Fig. 10 Biochar model predictions of corn grain yield, corn stover, soil bulk density (BD), soil pH, SOC (soil organic carbon) and vol- 
umetric soil moisture vs. experimental observations from Rogovska et al. (2014). The biochar was applied in November 2010 (see 
parameter values in Table 1) and the measurements were obtained in 2012. The simulation refers to a corn crop (Table 2) growing in 
a Clarion soil (Table 3) in Ames, IA. 


preliminary testing of the model. Comprehensive model 3. changes in CEC and pH that affect urea hydrolysis 
calibration and validation is the next step. Given the and nitrification rates (Eqns 11 and 12); 
diversity in biochar types and their complex interactions 4. changes in NH, pool due to mineralization/immobi- 
with crops, soils and climate, the task of calibrating and lization of N caused by the decomposition of biochar 
validating the biochar model must be iterative over time (Eqns 2, 7-9); 
and include multi-locations, multi-species and multi- 5. changes in NH, pool due to adsorption and desorp- 
biochar type trials. In that respect, the APSIM modeling tion caused by biochar CEC (Eqns 11, 13); 
platform becomes even more important as it is fre- 6. changes in field capacity and SW which affect water 
quently updated, freely available in the public domain stress factors and therefore indirectly affect all of the 
and contains several crop models and soil/environmen- N transformation processes (Fig. 4) and plant growth 
tal models along with flexible managers to assist investi- (Fig. 8); 
gation of biochar effects on diverse cropping systems in 7. changes in BD, which affect SAT and therefore 
different environments. WEPS. WEPS, in part, drives N20 emissions together 
The current version of the biochar model captures the with NO3, temperature and C in APSIM (Thorburn 
most important biochar effects on soils and extrapolates et al., 2010). 
the effects to crop yields by elaborating mechanistic Our next steps in biochar modeling are: (i) finalize 
crop models available in APSIM. Regarding the N on-going calibration tasks (e.g. Fig. 10) as well as 
cycling which is of particular importance, the biochar expand model calibration tasks to diverse environ- 
model affects N cycling via the following direct and ments, cropping systems and biochar types; (ii) 
indirect pathways (Fig. 1): develop a database with calibrated biochar parameters 
1. increase in the size of the SOC pools that can stimu- (Table 1) that will allow subsequent system analyses of 
late higher N mineralization rates (Eqn 2) biochar effects by the users; (iii) add new mechanisms 
2. changes in the turnover rates of FOM and SOC pools and improve existing ones in the model. For instance, 
due to priming (Eqns 4—6) quantification and modeling of biochar effects on root 
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functioning, seedling emergence, possible allelopathy 
effects, and distinction between fresh and aged biochar 
effects are some of the future challenges, among 
others. 


Knowledge gaps and research priorities to support 
mechanistic modeling 


In this study we encountered several knowledge gaps 
that would benefit from further research. We propose 
the following for prioritization: (i) NH, absorption and 
release patterns in the long term. The size of the NH, 
pool over time is the key to achieve production and envi- 
ronmental benefits (Fig. 1); (ii) experimental verification 
of the Saxton & Rawls (2006) pedotransfer functions for 
biochar response because correct prediction of soil 
hydrology affects the overall system substantially 
(Fig. 1); (iii) estimation of the C : N ratio of the labile 
and recalcitrant biochar separately because its value 
strongly affects the rate of N immobilization (Eqn 2); (iv) 
development of biochar specific temperature and SW 
modifiers; (v) experimental verification of the biochar C 
partitioning coefficients used in the model (Eqn 3; 
Table 1); (vi) quantification of the duration of priming, 
i.e. does the priming effect last for days, weeks, months, 
years and how is this related to biochar type, soil, and 
management?; (vii) distinction of fresh from aged bio- 
chars; (viii) the effects of fertilizers and residue manage- 
ment on pH (Yuan et al., 2011) needs to be addressed for 
a better representation of soil processes in the model. 

In summary, this study set the basis for a mechanistic 
and system assessment of biochar effects on soils and 
crops and their trade-offs. We envision that the biochar 
model will shed light on critical questions and generate 
hypotheses for experimental research. The model has a 
great potential to assist agronomists, soil scientists, 
economists, and policy makers to design sustainable 
and profitable biochar based cropping systems. 
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